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We calculate, using numerical methods, the Lyapounov exponent j(E) and the 

O , 

density of states p(E) at energy E of a one-dimensional non-Hermitian Schrodinger 

equation with off-diagonal disorder. For the particular case we consider, both j(E) 

> ■ 
ly-j . and p(E) depend only on the modulus of E. We find a pronounced maximum of 

p(|-E|) at energy E = 2/v3, which seems to be linked to the fixed point structure 

of an associated random map. We show how the density of states p(E) can be 

ov 

expanded in powers of E. We find p(\E\) = -4- + ^|i?| 2 + • • •. This expansion, 
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which seems to be asymptotic, can be carried out to an arbitrarily high order. 
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1 Introduction 

It has been realised over the last few years that a number of non-equilibrium problems can 
be described through the time evolution of a non-Hermitian random Hamiltonian. Some 
noteworthy applications include the motion of a particle in an imaginary vector potential 
U, vortex line pinning in superconductors ||, ||, and growth models in population biology 
||]. The interest for non-Hermitian random Schrodinger equations increased further 
U? §|, 0! Hi i> 0? Q) as it was reanse d that, at least in some cases, a transition from 



localised to delocalised states could take place even in one dimension. 

The simplest examples considered were one-dimensional tight-binding models for 
which the wave function satisfies 

e h ip n+1 + e"Vn-i + V n tfj n = Etfj nj (1.1) 

where V n is a random potential at site n and E is the energy. For real h and periodic 
boundary conditions (n = N + n), it was shown Jj], ||| that the eigenvalues are located 
(for a large system size) either on the real axis, or on lines which can be understood from 
the knowledge of the Lyapounov exponent for h = 0. 

Here we consider another simple one-dimensional non-Hermitian Schrodinger equa- 
tion, first introduced by Feinberg and Zee 0: 

e i9 "VWi+e ix "^-i=£^n. (1-2) 

In this case there is no random potential V n , and the disorder is purely off-diagonal. Both 
9 n and Xn are uniformly distributed between and 2-7T, and they are all independent. 

A numerical study of the spectrum of this Schrodinger equation suggested that the 
density of states was roughly uniform in a circle of radius ~ tt/2. We found it challenging 
to see whether existing tools JTB], |X3[ |T5], [HJ could be used to calculate this spectrum. 



In the present work we mostly study the density of states of Eq. (|1.2|) by calculating 
the associated Lyapounov exponent and using the Thouless formula [0 (which remains 
valid for non-Hermitian one-dimensional models as shown in Section |2|). 

The calculation of the Lyapounov exponent, and consequently the density of states, 
can be done either by Monte Carlo methods or by solving numerically an integral equation 
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for the Ricatti variable r n = \ip n /ip n -i\. This is done in Section |3|, where we also discuss 
the possibility of singularities both in the probability distribution of r n and in the density 
of eigenvalues E of Eq. (|1.2| ) . 

In Section [| we develop a method to perturbatively calculate the Lyapounov exponent 
and the density of states, in powers of the energy E. This approach is an extension of 
similar calculations done in the past, where the anomalous behaviour of the density of 
states at the band centre was obtained in the Hermitian case [ 15| , |I~6|1 . Our expansion 
indicates that the radius of convergence is E = and that the series is asymptotic. 



2 The Lyapounov exponent and the density of states 



An old and powerful way [13|, Q of studying one-dimensional random systems consists 
in rewriting recursion formulae like Eqs. ( |1.1| )-( [L~2"| ) in a matrix form. Indeed, Eq. ( |1.2j ) 
can be recast as 



where M„ is a random 2x2 matrix 



M n (E,9 n ,x, 



M„ 



4>n 
lpn-1 



(2-1) 



1 



(2.2) 



so that for any choice of the energy E and of the wave function ipo and ipi at sites and 
1, the wave function ip n for n > 2 reads 



4>n+l 



n 



Ee~ Wp —e^Xv-Op) 
1 



to 



(2.3) 



If ipo and %p\ are chosen arbitrarily {i.e. they are fixed, and not tuned functions of E and 
of all the 9 P and Xp)-> one expects that for large N 



(2.4) 



lim — log|^7v| =l(E), 

N—>co iV 



where j(E) is the largest Lyapounov exponent of the product of random matrices 

Because the matrices M n are invariant under the transformation M n (E, 6 n , x n ) ~ * 
M n (Ee l13 , 9 n + j3, Xn + P) for any j3, and because by shifting all the phases 6 n and Xn by 



a constant /3 one gets an equally likely random sample, it is clear that j(E) = ^(Ee^). 
Thus the Lyapounov exponent depends only on the modulus of E. For similar reasons, 
the average density p{E) of eigenvalues is invariant under the change E — > Ee % ^\ 

l{E)= 1 {\E\), 
p(E)=p(\E\). (2.5) 



In fact these two quantities are related by the Thouless formula 17 



/oo ^oo 

dE x dE y p(E x + iE y )log\E-(E x + iE y )\, (2.6) 

-co J— oo 

which can be derived as in the Hermitian case by the following argument: Consider a 
finite system of N sites with Dirichlet boundary conditions ("0o = ipN+i = 0). If one 
fixes ip = and ipi — 1, the ipw+i calculated from the recursion ( |2.1|) is found to be a 
polynomial of degree N in E, the zeroes of which are the N eigenvalues E a (for Dirichlet 
boundary conditions): 

/ N \ N 

ii> N+ i(E) = exp -i J2 On II ( E ~ E *)- ( 2 - 7 ) 

\ n =l / a=l 



Applying Eq. ( |2.4| ) for large iV leads to Eq. ( pi.6| ). 

It turns out that for all E ^ 0, the Lyapounov exponent j(E) is strictly positive 
by the Furstenberg theorem |jj|, [19| . (This will in fact be confirmed by our numerical 
results, and by the small- E expansion of Section ^.) One can then argue that the density 
of states is independent of the boundary conditions. Indeed, for general E, Dirichlet 
(ipN+i(E) = ipo(E) = 0) and periodic boundary conditions (ipN+i(E) = ipi(E)) are very 
similar: Both require that tpN+i(E) ~ e Nj ( E >, obtained by iterating Eq. (|2~3| ), takes 
very untypical small values at the eigenenergies. Thus, the eigenvalues corresponding 
to respectively periodic and Dirichlet boundary conditions should be very close, their 
distance being typically exponentially small in N. 

The Thouless formula ( |2.6[ ) can be inverted by using an analogy with two-dimensional 
electrostatics. If we interpret 'y(E) in Eq. (|2.6|) as the two-dimensional Coulomb potential 
created by a charge density p(E'), it follows from Poisson's equation that 

P (E, + IE,) = ± g + ^) TiE. + 1^), (2.8) 



or, exploiting the rotational symmetry of j(E) and of p(E), that 



P(\E\) 



+ 



2tt \d\E\ 2 \E\d\E 



l(\E\). 



(2.9) 



For \E\ > 2 the density of states vanishes. To see this, consider an eigenfunction ip n 
corresponding to an eigenvalue E a . By applying the triangular inequality to Eq. ( |1.2|) , 
one obtains \E a ip n \ < \ip n -i\ + iV'n+il- Choosing n such that \ip n \ = max (|^j|) i=1 , this 
implies \E a \ < 2. As the density of states p(E) vanishes for \E\ > 2, by expanding the 
logarithm in Eq. ( |2.6| ) and by using the rotational symmetry ( |2.5| ) of p{E) one readily 
finds that 

y(E) = log \E\ for \E\ > 2. (2.10) 



All our calculations which follow are based on another reformulation of Eq. (|1 . 2|) . 
Introducing the Ricatti variable r n 



A 



Ipn- 

the problem (|2.3| ) takes on the form of iterating a random map shown in Figure [1] 



(2.11) 



r n +i 



— + Ee 1 ^ 
r„ 



(2.12) 



where ip n (which is equal to Xn — 9 n + n modulo 2tt) is uniformly distributed between 
and 2ir. The Lyapounov exponent is then given by 



N 



l(E) 



N^oo /V , 
n=l 



(2.13) 



For large n the probability distribution of r n becomes independent of n, and the 
invariant distribution P(r, E) satisfies 

P(r, E) = r ^ r ds P(s, E)s(r- - + Ee ilp ) . (2.14) 

Jo 2tt Jo \ s ) 

From the knowledge of this invariant distribution, the calculation of the Lyapounov 
exponent j(E) follows from Eq. ( |2.13j ), and one has 

7(E) = / drP(r,£)logr, (2.15) 

Jo 

whereas the density of states p(E) is given by Eq. (|2[ 



All the rest of the paper is devoted to the determination of the invariant distribution 
P(r,E) solving Eq. Q2.14J) , and to the calculation of the Lyapounov exponent j(E) as 
well as the density of states p(E). 




Figure 1: Schematical illustration of the random map ( |2.12| ). Given r n , its iterate r n+ i is 
randomly distributed in the hatched region. The boundary constitutes two maps, T + and 
T_, obtained from Eq. ( |2.12| ) by choosing tp n = and tp n — ix respectively. For E < 2, 
these maps possess fixed points r + = T + (r + ) and r_ = T_(r^), that are respectively 
attractive and repulsive. For the particular energy E = 2/\/3 chosen on the figure, 
T(r i ) = r_. 



3 Numerical study 

In this section we numerically determine the invariant distribution P(r, E) solution of 
Eq. fl2.14| ), the Lyapounov exponent j(E), and the density of states p(E). 



According to Eq. ( |2.9| ) the computation of the density of states requires both 'j(E) 
and its first two derivatives. Whilst j(E) is rather easy to determine numerically, a 
precise estimate of its derivatives with respect to the energy turns out to be far more 
difficult. We have attacked the problem by two different methods. 

3.1 The Monte Carlo approach 

This first method consists in iterating Eq. ( |2.12| ) for a large sample. Typical shapes 
obtained for the stationary distribution P(r, E) are shown in Figure ^|. 

First, one notices that depending on the value of E, the support of the invariant 
measure P(r, E) is either finite of infinite. If one assumes that the support is an interval 
we see from Eq. (|2.12|) that if r n e [^min^max], the requirement that also 
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r n +i £ D"min, ''"max] gives conditions that r min and r max should satisfy. After analysing all 
the possible cases, one finds that the support of the distribution P(r, E) is the whole 
positive real axis for E < 2, whereas for E > 2 the support is finite, with r min and r max 
satisfying r min = E - l/r min and r max = E + l/r min , that is 

rmin = \ (E + VW^l) , r max = 1 (3E - VW^A) . (3.1) 

Second, for large enough E, there are apparent singularities at certain values of r. 
These singularities seem to be remarkable points of the maps T + and T_ obtained from 
Eq. ( |2.12[ ) by choosing tp n = or tp n = n, i.e. 

1 



T+(r) = E 

r 



Tjr) 



E- 1 - 

r 



(3.2) 



For each value r n , it is easy to see that T + (r n ) and T_{r n ) are the two extreme values 
that r n+1 can take: T_(r n ) < r n+l < T + (r n ). 
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Figure 2: Histograms of r n for energies E = 0.7, 1.2, 1.7 and 2.2, obtained using iV = 10 7 
iterations. Apparent singularities occur at the points r + and T_(r + ), marked by dashed 
lines on the figure. 
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Figure 3: Close-up of the most visible "singularity" of Figure |2].c, this time using N = 10 11 
iterations. 



The fixed point of T + 



E + VE 2 + 4 



(3.3) 



coincides with the most visible singularity of P(r, E) in Figure [| whereas the iterate 
T^(r + ) gives the position of another very visible singularity of P(r,E). Other apparent 
singularities seem to be located at the points T 2 (r + ) and T + (T_(r + )). 

In Appendix [A], we show that the singularities in P(r, E) are, if at all present, much 
weaker than they look on Figure |2|. This is confirmed by Figure ^|, which is an enlargement 
of Figure |2|.c, where there is a clear rounding of the singularity at r + . 

The calculation of the Lyapounov exponent is straightforward by Monte Carlo meth- 
ods. One just iterates the random recursion ( [2.12 ) a large number N of times (typically 
N = 10 8 ) and one obtains ^{E) by 



l{E) 



N+N' 



N 



E ioj 



(3.4) 



n=N> 



where N' is a large enough number of iterations to eliminate the transient effects due to 
the initial choice of r$ (here we took N' = N/10 which is exceedingly sufficient). 

For the density of states one needs Eq. (|2.9| ) to calculate the first two derivatives of 
■j(E) with respect to E. A possible approach would be to measure ^(E) from Eq. (|3.4| ) 
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for three nearby energies, E — AE, E, and E + AE, and calculate numerically the first 
and second derivatives. With this method it is however hard to find a good compromise 
for the choice of AE: If AE is too small, one does not have enough precision on the 
derivatives, and if AE is too large, all the structure of the density of states is artificially 
smoothed. 

To avoid this difficulty we tried to iterate directly the derivatives of r n with respect 
to E. The recursion (|2.12| ) has the form 

r n+1 = f (r n ,E,(p„). (3.5) 

Now, if r^ and r^ denote the first and the second derivatives of r n with respect to 
E, one can obtain (in a complicated form which we do not reproduce here because we 
wrote it in our programmes through a change of variables) recursion formulae for these 
derivatives 

r { nli = /i(r n ,r«, £,y? n ), (3.6) 

r£i = h{r n ,r£\r%\E,y n ). (3.7) 

Iterating these recursions N times yields an estimate for the first two derivatives oi^f(E) 



dE N *-L, TV 

n=N' 



2„jrn\ 1 N+N' 



d^(E) ^ 1 



- T 

AT ^-^ 



dE 2 N 

u~c iv n=N , 



r (2) /„(1)' 

n I n 



(3.8) 



Figure 3 shows our results for ■y(E) and p(E) obtained by this Monte Carlo approach. 
Whilst the results for j(E) seem to be alright, the density of states exhibits some irregu- 
larities. We tried to understand the origin of these irregularities by changing the length 
iV of the sample, the generator of random numbers, and the way in which Eqs. ( |3.6|) - 
(P-7|) were parametrised in our programmes. Although the positions of the irregularities 



were observed to change, we were unable to eliminate them altogether. Nevertheless we 
believe that they do have a purely numerical origin: From time to time, there is a very 
small value of r n produced by the iteration of Eqs. (|3.5|)-(|3~7|), and this gives such a huge 
contribution to the sums ( |3.8|) that all the remaining terms become negligible. Usually 
this big contribution is followed at the next step by another huge contribution of opposite 
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Figure 4: Lyapunov exponent ■y(E) and density of states p(E) obtained by Monte Carlo 
after 10 8 iterations. The inset in the left part of the figure shows the asymptotic approach 
to the exact result ( |2.10|) . In figure 4b, the vertical dashed line indicates the energy 
E = 2/V3. 

sign which more or less compensates the first big contribution. These events (where r n 
is very small) have a dramatic effect on the accuracy of our results, with the unfortunate 
consequence that the more we iterate, the more of these events we see, so that the more 
irregularities are visible. 

We unfortunately could not find a satisfactory way of avoiding these difficulties, and 
so we tried a completely different approach discussed in the next subsection. 

In Figure f|.b, we see that the density of states p(E) exhibits a non-trivial structure 
with a pronounced maximum at an energy E ~ 1.15. Trying to understand the origin of 
this maximum, we realised that at energy 

2 



E = —= ~ 1.1547- 
\/3 



(3.9) 



the point T_(r + ) becomes a fixed point of TL. In a more physical language, if one 
considers an infinite sample for which 



6 n = and Xn = ^ for n < 0, 
6 n = and Xn = for n = 0, 



6 n = Ti and Xn = for n > 0, 



(3.10) 
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one finds that there is a bound state at energy E = 2/y / 3 corresponding to a wave 
function ip n oc 3~^ n ^ 2 in addition to the continuous part of the spectrum at imaginary 
energies, E = ix with —2 < x < 2. 

We believe that the origin of the sharp maximum in p(E) at energy E = 2/v3 is due 
to the existence of these bound states, as in the case of impurity bands for Hermitian 
problems fl9|, |20fl . We did not succeed, however, to determine whether the density of 
states is analytic at this special energy E = 2/\/3, or whether it presents some weak 
non-analyticity. 

Our numerical results for ^{E) in Figure f|.a can be fitted for small E by an expression 
of the form 

7 (£) = k x E 2 + k 2 E* + ... (3.11) 

with k\ = 0.1595 ± 0.0005 and k 2 = 0.017 ± 0.001. To be more precise, we first find 
ki as the extrapolated interception of ^(E)/ E 2 with the ordinate. Then, subtracting 
the k\E 2 term, k 2 can be similarly found from the residual E A dependence. Needless 
to say, this procedure greatly enhances numerical noise at each stage, and so we had 
to content ourselves finding the first two terms. These are in good agreement with the 
results k\ = -^ ~ 0.1592 and k 2 = ^ — 0.01689 of the perturbation theory to be 
developed in Section f| (see Eq. ( |4.25| )). 

3.2 Discretisation of the integral equation for P(r, E) 

To avoid the problems of the Monte Carlo method, we tried to numerically solve the 
integral equation Q2.14Q for P(r, E). To this end the r- variable was discretised by taking 



1000 values of r at the points r k = 2fc/(2001 - 2k) for 1 < k < 1000. The integral 
operator ( |2.14| ) then becomes a 1000 x 1000 matrix, and the stationary distribution can 
be obtained by simply iterating a linear system. 

Figure § shows distributions P(r, E) obtained that way. They are very similar to what 
was obtained by the MonteCarlo method, cfr. Figure @.b-c. In Figure |6] we display the 
corresponding Lyapounov exponent 'y(E) and the density of states p{E). The Lyapounov 
exponent is obtained simply by the discretised version of Eq. (|2.15|) , and its derivatives 
are then calculated by replacing them by finite differences with a differentiation interval 
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Figure 5: Invariant distributions P(r,E) for energies E = 1.2 and 1.7, obtained by 
discretisation of the integral equation ( |2.14| ). 

AE = 0.025. (This value results from a compromise: For smaller AE a complicated 
structure appears in the derivatives due to the discretisation of r, whereas a too large 
AE smears out details of p{E).) 

We see that both 'j(E) and p{E) have essentially the same shapes as we have previ- 
ously encountered in figure 4, using the Monte Carlo method. The irregularities which 
were visible in Figure f|.b have now disappeared in Figure |6|.b. Moreover, the inset in 
Figure |^. a, showing the asymptotic approach towards the exact result ( |2.10[ ), indicates 
that the accuracy of ^(E) has been improved by roughly two orders of magnitude. 

In the insets of Figures |Ja and §.a, the difference j(E) — log(E) seems to vanish at 
E ~ 1.8. We believe that this difference is non-zero up to E = 2 but that it becomes 
very small (p(E) has an essential singularity: it vanishes and all its derivatives vanish 
at E = 2). The argument that eigenstates exist up to E = 2 can be borrowed from the 
theory of Lifshitz tails [2C, |IH| in the Hermitian case: in a random sample, one can find 
arbitrarily large regions where 9 n ~ Xn — 0- 



4 Perturbation theory 

In this section we develop a perturbation theory in powers of E to determine the in- 
variant distribution P(r,E) solution of Eq. (|2.14j ). When this is known, we can obtain 
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Figure 6: Lyapunov exponent j(E) and density of states p{E) obtained by discretisation. 
The irregularities visible in Figure f|.b have now disappeared. 



the corresponding perturbative series for 'y(E) from Eq. (|2.15|) , and that of p{E) from 
Eq. O. 

For E = there are infinitely many solutions to Eq. (|2.14|) . Namely, the iteration 
formula fl2.12| ) reduces to 

r„+i = -, (4.1) 



and any distribution preverving the symmetry r <^ 1/r (i.e. such that P(r, 0) = P(l/r, 0)/r 2 ) 
is a solution. From Eq. (|2.15|) we thus immediately have 



T (E = 0) = 0. 



(4.2) 



On the other hand, for \E\ > the invariant distribution P(r,E) is unique, as wit- 
nessed by the numerical results of Section [j[ So the expansion in powers of E that we 
are going to develop presently is a degenerate perturbation theory, as there are many 
solutions when E = and a single one, P(r, E), when E > 0. 

It is useful to analyse first the small-i? expansion of the integral operator which 
appears in Eq. ( |2.14|) . Consider two distributions Q(r) and Q(r) related by 



Q(r)= C^ rdsQ(s)5[r- 

JO Z7T JO 



- + Ee lv 
s 



(4.3) 



13 



When E < r, this can be rewritten (by integrating over s) as 



Q(r) 



2w dip r d + rE 2 — 2rE 2 sin if + IrE cos if J r 2 — E 2 sin if IE cos <f + Jr 2 — E 2 sin if 



2vr 



E 2 ) 2 J r 2 -E 2 sin 2 if 



Q 



-E 2 

(4.4) 



Assuming E <C r the integrand can then be expanded as a power series in £?, and after 
performing the integral over if we find that 



Q{r) 



;Q 



+ 



E 

T 



2 r 



9 



Q 



^' J 



M; 



64 



225 



+ 



E 



6 



11025 



351 



\ 1" / 1"' I T 



^e" f i 



22 



2304 
732 



7 "' 



,,.1\ 25839 _,/l 

<2 - + — 5-Q' 



,12 







(4) 



+ 



45 

7l3 



r 



9 



18261 A, 



iy _LU \ r f } 1™ 



+ 



1 

5382 



5^ (3) ~ +^ (4) '~ 



-g 



(3) 



+ 







(5) 



U +^ (6) 



0(E 8 ). 



(4.5) 



We see that the effect of a non-zero energy is that Q(r) is not just trivially related 
to Q (M, but also depends on all the derivatives of Q- Quite remarkably, the above 
expansion can be written in a compact form valid to all orders. Defining the second 
order differential operator M by 



the expansion of Eq. 



*, d d 1 
dr dr r 

to an arbirary order in E can be written 



oo ZT^n 

Q(r) = £ —^M r 



n=0 



n! 



V 



(4.6) 



(4.7) 



It is straightforward to see that this agrees with the naive expansion formula ( |4.5|) . The 
validity of Eq. ( |4.7| ) to an arbitrary order in E is established in Appendix B. Note that 
a simple change of variable r — > 1 jr transforms Eq. ( f4.7| ) into 



Q 



where the operator L is defined by 



E 

n=0 



E 



2n 



A n (n!) 



:L n 



r 2 Q(r) 



(4- 



d d 



L = r z — r — r. (4.9) 

dr dr 

Coming back to the invariant distribution P(r,E), we see by combining Eqs. 



and ( |4.8j ) that it satisfies 



oo / rp2 \ n n -i 

P(r ' E) = SW p ?o[P'(n-p)!] 2 



M p 



-L 



71— P 



(r 2 P(r, £) 



(4.10) 
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If we look for a solution P(r,E) which can be expanded as 



E 



/£2\2 



P(r,E) = P (r) + —P 1 (r)+l — \ P 2 (r) 



(4-11) 



we obtain by equating the two sides of Eq. ( |4.1U| ), order by order in P, that 

1 
M v 

— ' In! In — TillM 
n=0 p= 



m+l n | 

Pm+l(r) ^?PF* 



■ n— p 



( r "m+l-n 



_*, r 



(4.12) 



The term with n = on the right-hand side is just P m+ i(r), same as the left-hand side, 
so that P m+ i(r) disappears from the equation. Therefore 



m+l n -i 

n ?i p ?ob!(«-p)!] 2 



M p 



L""P (r 2 P m+1 _ n (? 



0. 



(4.13) 



Moving all but the two n = 1 terms, which are the only ones involving P m (r), to the 
right-hand side we arrive at the equation 



M 



:i+np m (r] 



m+l n 
n=2 p=l 


1 

[p! (n — 


pW 


MP- 1 


m+l i 


r 2rn-l 


(r 2 P 


.,1 (r 



;L n ' p (r 2 P n 



+i-n(r 



(4.14) 



Note that an M operator has been factored out (by using the fact that Lr 2 = r 2 Mr 4 ) 
from all the terms. Integrating over this, and shifting the summation variables, we finally 
find 



l + r 4 )P m (r) 



1 



^rtogr + B m r-E- [(n+1) , ]2 



EE 



i 



ti^ol(P+^(n-p)\? 
where A m and B m are a priori arbitrary constants. 
When m = 0we find 

r logr 



M p 



-n—p 



L \T P m -n\> 

(r 2 P m ~n(r) 



(4.15) 



Pn r) = Ac 



l + r l 



+ B 



1 + r 4 



(4.16) 



By inserting this into Eq. ( f4.7| ), with Q(r) and Q(r) replaced by P(r, P), we see that for 
Po(r) to satisfy this equation to zeroth order in E, i.e. Po(r) = P (l/r)/r 2 , the constant 
A should vanish, whereas B could be arbitrary (the problem is linear and B would be 
fixed by the normalisation of P(r, E)). This leads to Po(r) = B r/(1 + r 4 ). 
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For the same reason, namely that P(r } E) is the fixed point of Eq. (O) rather than 
the solution of Eq. ( |4.1(J| ) obtained by by iterating Eq. ( |4.7| ) twice, one can show that 
A m = and B m is arbitrary for all m. In fact, a non-zero B m produces a contribution to 
P m (r) proportional to Po(r), and so one might as well choose B m = for all m ^ by 
allowing B to depend on E. In practise, then, all the A m = for m > 0, all the B m = 
for m > 0, and B is an arbitrary function of E (which can be fixed by requiring that 
the distribution P(r,E) be normalised). 

In this way one arrives at 

-.3 Af^ 



P(r,E) = Bo(E) 



l + r 4+E ^(1 + r 4 ) 2 (l + r 4 ) 4 J +C> ^ 



(4.17) 



The recursion ( ^4.15|) gives an explicit formula for the invariant distribution to all 



orders in powers of E. We now show how to produce the equivalent expansion for the 
Lyapunov exponent. First, it is easy to check by induction that all the P m (r) are finite 
sums of functions of the form 



,3 



fn(f) = -, r^ and q n (r) = -— . (4.18) 

Jny - I n _i_ r i\n i " lv I (1 4- r) n 

Indeed, this is true for m — 0, and higher orders are generated from Eq. (|4.15|) by means 
of the operators M and L. These operators act on f n and g n as follows: 

Mf n = I6n 2 g n+1 - 16n(n + l)g n +2, 

Mg n = 4{l-2n) 2 f n -32n 2 f n+1 + 16n{n + l)f n+2 , 

Lf n = A(l-2n) 2 g n -32n 2 g n+1 + 16n{n + l)g n+2 , (4.19) 

Lg n = 16(1 -n) 2 /„-i- 16(1 -3n + 3ri 2 )/ n + 48n 2 /n+i-16n(n + l)/ n+2 . 

Applying Eq. ( f4.15| ) repeatedly, we have found the explicit expressions for the first seven 
orders of the (unnormalised) invariant distribution. We report the first few orders here, 
and defer the remaining ones to Appendix y. 

Po = fi, 

-j = 92-4g 4 , 

p i 30 

-^ = 2/2 + y h ~ 22/4 - 158/ 5 + 320/ 6 - 160/ 7 , 
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P 3 1 673 131 4456 842 MoAn 

64 = 2 92 ~ l8 93 ~ l8 9i + — 95 ~ ^ 9Q - 2lM ° 97 

+ 56520^ 8 -53760^ 9 + 17920^io. (4.20) 

In order to compute the contributions to j(E) and to the normalisation, one needs the 
following integrals 



o 

oo 



/0 

Using ( |4.18| ), these are easily evaluated 

2n - 1 T . . T 7i 

*n+l = ~^Z ■*« Wlth -'I = T> 



dr/ n (r), 

/•oo 

J n = / drg n (r), 

Jo 

/•oo 

K n = / dr/ n (r)logr, 
Jo 

POO 

L n = / dr^ n (r)logr. (4.21) 



2n 


-'n 


1 




4(n - 


-1)' 


2n- 


1 K 


2n 


J^-n 


71—1 


T, 



4(n — i) 

K n+1 = ^— — -K n - —I n with K x = 0, 

An 

L n+1 = - — -L n - — J n with L 2 = 0. (4.22) 

n An 



Combining all of this we arrive at the expansion 



^ = ,G(EHH(E) < 4 ' 23 ' 



with 



F(E) = ±E 2 + —E e + —E l0 -^-E 1A + O(E ls ), 
K J 8 288 3600 235200 v ; 

G{E) = 1 + _*_& + ^gL g8 + «^«l ^ + 0(£W) , (4.24) 

v ; 4 256 294912 56623104 v ; v y 

ffW = _i.fi. _ ___.fi. + 1088W _■■■ + 576717747329 

v ; 12 15120 29937600 490377888000 v ; 

It is worth noting that all the coefficients in the expansions of the functions F, G and H 
are rational. Equivalent ly, 

liE) _ ±& + _L E < + ( i ____W + ( r i — *L.)_« 

' V ^ 27T 67T 2 V187T 3 11527T/ V547T 4 201607T 2 / 
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Figure 7: Log- log plot of |o2n| 1//< - 2n * ) versus 1/n. 
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(4.25) 



In this way we can in principle generate j(E) to any order. We calculated numerically 
the coefficients a2 n of the expansion of ~f(E) 

7 (E) = a 2 E 2 + a^E 4 + ■■■ + a 2n E 2n + ■■■ 

up to n = 46 in order to estimate the radius of convergence of this series. In Figure |7| 
we show a log-log plot of |a2n| _1 ^ 2n ^ versus 1/n. For n > 10 the data are well fitted 
by the form \a 2n \~ 1 ^ 2n ^ oc rr v with v ~ 1.08, shown as a dashed line on the figure. 
This result indicates rather convincingly that the large-n limit is zero, and thus that the 
radius of convergence vanishes. We therefore believe that our expansion is an asymptotic 
expansion. 



18 



5 Discussion 

In this paper we have determined numerically the Lyapounov exponent j(E) and the 



density of states p(E) of the one- dimensional non-Hermitian Schrodinger equation (1.2) 
when the phases 8 n and \ n are uniformly distributed. We have also developped a method 
of expanding these quantities in powers of the energy E and obtained j(E) up to order 
E 16 . Our procedure can be extended to obtain, in principle, an arbitrarily high order. 

The expansion in powers of E is based on the expansion of the integral equation 
( p,14j ) satisfied by the invariant distribution P(r, E). Since this expansion is in principle 
only valid for E <C r, there is a possibility that the perturbation series for j(E) suffers 
from the fact that when we integrate over r, we use an expression valid for r ^> E even 
in the neighborhood of r = 0. One cannot exclude for example that when r and E 
become comparable (and small), P(r,E) becomes a scaling function of r/E. This could 
invalidate our expoansion of ~f(E). However when we compared our expansion with the 
results of the simulation we found an excellent agreement for the first two terms and so 
it is possible that our expansion, a priori valid for r 3> E, remains valid down to r = 0. 

Our numerical results show that the density of states vanishes outside a circle of 
radius \E\ = 2, is non-zero inside this circle \E\ < 2, and has a non-trivial shape with 
a pronounced maximum at an energy E = 2/a/3. The analyticity of j(E) or of p(E) at 
this energy remains an open problem. Another interesting question would be to predict 
how p(E) vanishes at E = 2. 

One could try to extend our approach to other situations, in particular to cases where 
the phases 8 n and \ n have a non-uniform distribution ||. In that case, the phase and 
the modulus of the ratio ip n+ i/ip n are in general correlated and one should look for the 
invariant distribution P(r, E) of the complex variables r n and E. Still, the fact that 
^{E = 0) = would remain true, and one could try to expand the invariant distribution 
in powers of E and use the Thouless formula to calculate p(E). 
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A Nature of the singularities of P(r 1 E) 

In this Appendix we discuss the occurence of singularities in the invariant measure P(r, E) 
at the fixed point r + of the map T + (r) = E + -, when one iterates the random map ( |2.12|) 

' ' : " (A.1) 



+ Ee 1Lpr 



r n +i 

with uniformly distributed tp n between and 2n. 

We believe that the apparent singularities in P(r, E) seen in Figure || are due to 
events where several successive r n are close to the fixed point r + , and that the nature of 
the singularity can be understood by analysing only the neighbourhood of r + . 

In physical terms, there is a competition between the deterministic part of the map 
T + (r) which tries to concentrate the distribution on its attractive fixed point r + , and the 
noise due to ip n which tends to broaden the distribution. To analyse the neighbourhood 
of r + , we consider the linearised problem 

s n+l = ~ a S n + t n , (A-Zj 

where a e]0, 1[ is a fixed slope, and t n is a random positive number. The variable s n 
represents the difference r + — r n , when this difference is small. By expanding Eq. ( |A.1| ) 
to second order in <p n , it is seen that for Eqs. ( |A.1| ) and ( |A.2|) to be equivalent in the 
neighbourhood of r+, one should choose 

g + 2-WgTI 

a = -T + (r+) = , 

E 



"» 2(Er + + 1) ^' (A ' 3) 

The essential feature of the distribution of t n is that all the t n are positive and that 
the distribution Q(t) of t n diverges as Q(t) ~ t -1 ^ 2 as t — * 0. If we iterate ( |A.2j ) 
numerically for E = 1.7, that is for a = 0.213851 • • •, we recover at s = a singularity (see 
Figure |8].a) which ressembles the one seen on Figure ^|.c. Actually, the two distributions 
are approximately mirror images near their respective fixed points, but this is due to the 
fact that s n ~ r + — r n rather than s n ~ r n — r + . 

Trying to analyse the nature of this singularity, we see from Eq. ( |A.2| ) that s n can be 
written as 

s n = t„-i - a£ n _ 2 + a 2 t n -3 - a 3 £«-4 + aH n - 5 - . . . (A.4) 
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Figure 8: Invariant distribution of s n for the linearised problem defined by Eqs. ( |A.2| )- 
( |A.3| ), here with E = 1.7. The magnification (obtained with N = 10 11 iterations) shows 
that there is in fact no singularity at s = 0. 

or as 

s n = y n - ay' n , (A. 5) 

where y n and y' n are two positive independent and equally distributed random variables 
of the form 



Vr, 



T n + a r n _i + a t„_ 2 + a r„_ 3 + • • • , 



y'n = < + « 2 Cl+« 4 <-2 + « 6 <-3 + ---- (A.6) 

All the r n and the r' n are independent and distributed according to the same distribution 
Q(t) as for the t n . 

Let us first consider the case where s n is still given by Eq. ( |A.5| ) but where the 
distribution ir(y) of the positive variables y is arbitrary. Using Eq. ( |A.5| ) one can show 
that the distribution P(s) of s is given by 

("OO 

P(s) = / dyrr(ay + s) n(y) for s > 0, (A. 7) 

Jo 

P(s) = / dy7c(ay)7c (y - -j for s < 0. (A. 8) 

One can calculate the successive derivatives with respect to s at s = 0, and one obtains 
for the n th derivative 



p(«)(0 H 



dyir {n) (ay)7r(y), 



(A.9) 
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p(»)(0_) 



dy7i(ay)7i in \y). 



(A.10) 



If we assume that ir(y) and its first n ^- derivatives vanish at the origin, and as y — > oo, 
it is easy to see by integration by parts that Eqs. (|A.9|) and (|A.10|) coincide, so that P(s) 
has at least n derivatives at 0. 



We are now going to argue that since y n is given by Eq. ( A.6 ), its stationary distri- 
bution ir(y) vanishes at y = as well as all its derivatives. To see this, let us consider 
the generating function (e _/3y ) of y. From Eq. ( |A.6| ) one has 



°Q r roo 

- ( *) = U / dtQ(t)e 

n=0 LJ0 



-f3ta 2 



(A.ll) 



Because Q(t) ~ t 1 ^ 2 for small t, one has for large (3 

For large (3, one has (e'^) ~ (e _/3a y ) and one finds that to leading order 

llog 2 /5 



(A.12) 



log (e"«> 



8 log a 



(A.13) 



We see that (e /3i/ ) vanishes faster (as log a < 0) than any negative power of (3 as f3 — > oo. 
As a consequence ir(y) and all its derivatives vanish at y = 0. 

Therefore we can conclude that all the derivatives of the distribution P(s) are con- 
tinuous at s = 0. This is confirmed by a magnification of the small-s region, as seen in 
Figure |8].b, where the rounding of P(s) at s = becomes visible. 

Coming back to the stationary distribution P(r, E), we observe the same rounding of 
the apparent singularity at r + (see Figure |3|) . 

B Justification of Eq. ( |3t71 ) 

Imagine that s is a random positive variable distributed according to some given distri- 
bution Q(s) and let r be given by 

1 



+ Ee lip 



B.1 
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where <p is uniformly distributed between and 2n. We want to show that the distribution 
Q(r) of r can be written as in Eq. ( |4.7| ): 



oo pf'tn r 1 / 1 \ 



(B.2) 



n=o 4" (n\y 

where 

M=^r^i (B.3) 

dr dr r 

Let us assume for simplicity that all (positive and negative) moments of Q(s) exist. 
By taking the (2p) th power of Eq. ( |B.1| ) and by averaging over tp, one can see that 

We see that for Eqs. ( |B.2| ) and ( |B.4| ) to be equivalent, one simply need to show that 

W Tds — 

[{p-n)\} 2 J s 2 p~ 

For n = this is an obvious identity, and for n > it can be established by induction, 
using integrations by parts. 



C Higher orders of P(r, E) 

Here we list the higher order terms of P(r,E), obtained from Eq. (|4.15| ), expressed in 
terms of the elementary functions f n and g n . These higher order terms are used together 
with Eq. ( 4.20J ) to obtain the expansion ( 4.23 ) for the Lyapunov exponents. 



1 r°° r 1 /1\i Ml 2 roc i 

- / drr 2p M n -Q - = .. m — / ds — — Q(s). (B.5) 

\: n Jo [r 2 \rj\ \(p-n)\] 2 Jo s 2 P -2n^^ > v > 



P 4 221 16295, 94193 2529101 2619745 

44 - ^/2 + ^-/ 3 -^-/4-25803/ 5 + h + g h 

_ 22135892 ^ + 77875796 ^ _ 43988480/io + 40636800/ii _ i 97 i2000/ 12 

+ 3942400/ 13 , (C.l) 

P 5 413 370801 11887541 34137437 76908091 

45 72 yi 64g va 324Q y 4 135 y & 13g yb 

174769499 684150977 1394447776 38774321284 

-97 H o ^8 H o ^9 « S'io 



9 9 3 9 

+ 14632O268160H - 27725662640^i 2 + 31886412800^i 3 - 22180787200^14 

+ 8610201600^15-1435033600^16, (C.2) 
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P 6 _ 132619 55698437 _ 1031521993 10179687761 76720899118 

4^ ~ 288 ^ + 2592 * 3 1080 ^ 1080 + 405 ^ 

166044285683 2068679706428 3096140849561 1697546657546 

+ 270 h 135 /8+ 135 h+ 3 h ° 

147448995026264 „ 

- 4380421880262/n + f 12 - 38217908241192/i 3 

9 

179597752851712 192990668440960 , . avmno .. fl1 . n .. 
+ /i4 /i5 + 46970084761600/ 16 

22329532825600/ir + 6245266227200/i 8 - 780658278400/i 9 (C.3) 

P 7 _ 288203 1354843111 815716196107 16419260711323 
47 ~ 2592 92 23328 9s + 583200 9i + 170100 ^ 5 

207216071509033 113911389465407 579469587812377 



9% z~—z 9i H z~—z 9s 



340200 u 5670 "" 5670 

677054300749364 1078075200453799 324458646218758 
+ Wb 99 ri 9io g <7n 

8114117353035130 50273274830679712 180446654936149976 

912 ~ #13 H n ^14 



9 " 9 " 9 

433731592050615872 735643743703812320 298906446014525440 

9i5 H 7, 9w = 9n 



9 J 9 J 3 

+ 2612301 6 31 11200 ^ 18 _ 53361562052198400^ 9 + 21803255983308800# 20 
o 

- 5339702624256000^ 2 i + 593300291584000# 22 . (C.4) 
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